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We have carried out a coordinated experimental and the- 
oretical study of single-electron traps based on submicron 
metallic (aluminum) islands and Al/AlO^/Al tunnel junc- 
tions. The results of geometrical modeling using a modified 
version of MIT's FastCap were used as input data for the 
general-purpose single-electron circuit simulator MOSES. The 
analysis indicates reasonable quantitative agreement between 
theory and experiment for those trap characteristics which 
are not affected by random offset charges. The observed dif- 
ferences (ranging from a few to fifty percent) can be readily 
explained by the uncertainty in the exact geometry of the 
experimental nanostructures. 
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I. INTRODUCTION 

Recent advances in the physics of single-electron charg- 
ing of macroscopic conductors (for general reviews see, 
e.g., Refs. HID) have led to proposals for several new 
analog and digital electronic devices. Such devices are 
considered, in particular, to be the most likely candi- 
dates to replace silicon transistors in future ultra-dense 
electronic circuits - see, e.g., Refs. HQ]. 

Single-electronics is presently one of the most active 
areas of solid state physics and electronics, with hundreds 
of experimental and theoretical works being published 
annually. We are not aware, however, of any previous 
attempts to quantitatively compare experimental data 
for a particular device with results of theoretical analysis 
including geometrical modeling 0|. Such a comparison 
was the main objective of this work. To that end, we 
selected one of the simplest devices, the single-electron 
trap|§ 

Figure u\ shows the schematic layout of the circuit we 
discuss in this paper, which consists of a trap coupled to 
a single-electron electrometer. We will distinguish two 
types of conductors ("nodes") in the circuit: externals, 
wires which extend to the edges of the chip and con- 
nect to the external measuring devices; and islands, small 
metallic segments that are connected each other and to 
the externals by tunnel junctions. 
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FIG. 1. Schematic view of the 8-junction single-electron 
trap/electrometer circuit. Islands 12 and 13 are strongly cou- 
pled and together provide the energy well for an extra elec- 
tron. Islands 6-11 form the array separating the well from the 
drive external 1. 

The trap consists of a larger island, providing the 
potential well for the extra electron, separated from a 
voltage-biased "drive" external by an A^island array. 
The islands of the array are linked by (N — 1) tunnel 
junctions with low capacitance Cj and tunnel conduc- 
tance G,-: 
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Under condition (g), each electron is localized inside a 
single island at any given time. As Fig. ||a shows, the ar- 
ray creates an electrostatic energy barrier AW ~ e 2 /Cj 
between the drive electrode and the trap island. To in- 
ject an additional electron into the trap, a bias voltage 

V = Vi — Vi is applied to the device. At a certain value 

V = y+ the energy barrier is suppressed: an electron 
tunnels from the drive external through the array and 
into the trap island. To extract the electron from the 
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trap island, a voltage V = V- is applied, causing a hole 
to tunnel from the drive external to the trap island, an- 
nihilating the trapped electron. 
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FIG. 2. Energy of a typical circuit calculated for three val- 
ues of Vi (V+, V eq , and V_), with V2 = 8 mV: (a) as a function 
of the position of one electron in the array; (b) as a function 
of the position of one hole in the array, with one electron in 
the trap (node 13). 

Electrons can also overcome the energy barrier by ther- 
mal activation and by macroscopic quantum tunneling of 
charge ("cotunneling"). At sufficiently low temperatures 
(|l|) , the rate of thermally activated hopping over the bar- 
rier is roughly 
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while the rate of spontaneous cotunneling through the 
barrier scales as 
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If conditions ([!]) and (||) are satisfied, and the number 
N of junctions in the array is large enough, the rates 
of thermal activation and cotunneling may be very low. 
Thus, the lifetime 77, = (Ft + Tq) 1 of both the zero- 
electron and the one-electron states of the trap may be 
quite long, and the device may be considered bistable. 

When the voltage V is driven beyond the threshold 
V+ or V-, the electron or hole tunnels through the ar- 
ray in time r ~ Cj/Gj, which may be many orders of 
magnitude shorter than r^. Thus, in principle, the trap 
can serve as a memory cell. Its contents can be read out 
non-destructively by capacitive coupling of the trap to 
the single-electron electrometer |l]-|^] (see Section 0). 

Early attempts to trap single electrons were made 
by Fulton et al. ||, using systems with two and four 
Al/AlOj; junctions of area ~ 100 x 100 nm 2 at tem- 
peratures down to 0.3 K. Their results implied trapping 
times tl ~ 1 sec. Similar experiments by Lafarge et 
al. H yielded tl < 1 sec, much shorter than could be 
anticipated from formulas (||) and (||). A later attempt 
Q used a semiconductor (GaAs) structure with a nar- 
row 2DEG channel instead of a well-defined tunnel junc- 
tion array. A bistability loop was observed, but its size 
was not clearly quantized, implying that the number of 
trapped electrons was much larger than one (the authors 
estimated this number to be 80-100). 

Finally, Al/AlO^ trap circuits designed and fabricated 
at Stony Brook |l(],[ll| yielded trapping times of over 10 4 
sec (limited only by observation time). The main goal of 
the present work was to compare the experimental data 
obtained for these traps with a quantitative theoretical 
analysis of the circuits. For this purpose, we have con- 
structed a geometrical model of the circuit, calculated 
the full matrix of self- and mutual capacitances for the 
conducting nodes in the model, and simulated static and 
dynamic properties of the trap using these capacitances. 



II. FABRICATION 

Circuits consisting of two layers of partially overlap- 
ping nodes were fabricated using the standard shadow 
mask technique fl^jis]] . The process begins with a Si 
substrate, either stripped of oxide or covered by a layer 
of Si02 of thickness H=500 nm. The substrate is coated 
with a PMMA/copolymer double layer mask. The cir- 
cuit pattern is written onto the mask using a scanning 
electron microscope. Then the mask is developed, the 
Al circuit elements are deposited onto the substrate, and 
the mask is lifted off. The fabrication process is described 
fully in Ref. KW. Here we present the essential details. 



A. Mask 

The circuit layout, consisting of a set of line segments, 
is first specified in a "mask file" (Fig.|). A version of the 
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same file is also used to start the computational modeling 
process (see Section [II). Wide lines in Fig. || represent 
the parts of the externals that extend from the trap and 
electrometer to contact pads at the edge of the chip. The 
narrow lines extending inward from the wide lines (Fig. ^) 
represent the inner parts of the externals. The short, 
narrow line segments (Fig. ^a) represent islands. 
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FIG. 3. Mask for complete chip containing several circuits. 
Pattern for circuits discussed in this paper is circled. 





FIG. 5. (a) Central part of mask, (b) AFM image of central 
part of fabricated circuit, (c) 3-D outline of central part of 
geometrical model used for capacitance calculation. 

The pattern of lines is written on the mask using the 
electron beam. Upon chemical development, each line 
in the PMMA becomes a window opening into a larger 
cavity in the copolymer, which is more susceptible to 
the electrons. This procedure results in a mask, shown 
schematically in Fig. [j], with several suspended bridges. 



FIG. 4. Closeup of layout of circuits discussed in this paper, 
showing 20 /J,m cutoff radius used in simulations, and external 
node numbers. Externals have wide (W = lfim) and narrow 
(w ~ 50 nm) sections. 
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FIG. 6. Schematic view of shadow mask fabrication: (a) 
first-layer evaporation, (b) oxidation (oxide thickness exag- 
gerated for illustration purposes), (c) second layer evapora- 
tion. 



B. Deposition 

The aluminum islands are deposited in two layers. The 
first layer is resistively evaporated in high vacuum di- 
rectly onto the room-temperature substrate (Fig. |6|a). 
This layer is then oxidized at ~ 10 mTorr O2 for ~ 10 
min (Fig. ^b), covering the Al islands with a ~ 1 nm 
layer of A10 x . Before depositing the second layer, the 
chamber is re-evacuated and the substrate is tilted rela- 
tive to the Al source. The tilt creates a shift s between 
these two groups of islands, so that they partially over- 
lap (Fig. ^c). The A10 x creates tunnel barriers between 
the first and second layer islands. In our circuits, the 
shift s was about 120 nm along the vertical direction in 
Fig. 0a, b. The second aluminum layer is made thicker 
than the first, to allow reliable step coverage. 

An AFM image of the resulting circuit is shown in 
Fig. 0b. This image exaggerates the island widths be- 
cause of the finite angle of the AFM tip. Other obser- 
vations (including SEM imaging) show that the islands 
oriented perpendicular to the direction of the shift were 
in fact spatially separated, in the successful samples. 

Figure 0c shows a simplified model of the central 
part of the circuit, with externals and islands numbered. 
There are two islands for each corresponding window in 
the mask. For example, islands 6 and 7 are the first- 
and second-layer products of the same window (see also 
Fig. 0c.). Since the two layers of each external overlap 



each other extensively and are connected to the same 
voltage/current source, they effectively serve as one con- 
ductor. Thus, there is only one external for each corre- 
sponding window in the mask. 



III. GEOMETRICAL MODELING 

The essential electrostatics of a group of conductors 
can be described by their mutual capacitance matrix, C. 
A program known as FastCap |Q can calculate C for an 
arbitrary collection of conductors, given the geometry of 
the conductors as input. The conductor surfaces are pre- 
sented to FastCap as a set of discrete elements, or "pan- 
els" . We wrote a program called Conpan (for conductor 
panels) to generate a 3D paneling of a simplified model 
of the experimental system, starting from a 2D mask file. 
We will first explain the Conpan algorithm, then how its 
input parameters were derived from experiments. 



A. Conpan Algorithm 

Conpan represents circuit nodes by means of data 
structures called "sections" . Each section is a collection 
of data about a node or part of a node. The data include 
parameters such as node number, layer number, and lim- 
its in the xy plane. Sections may be recursively divided 
into subsections to represent overlaps and to facilitate 
paneling. 

Consider two line segments from the larger mask file 
(Fig. ^a). These two segments eventually produce four 
islands separated by three tunnel junctions. Conpan ex- 
pands each segment into a first-layer section (Fig. 0b) us- 
ing the line-width w. The second-layer sections (Fig. 0c) 
are initially identical to the first-layer sections except for 
a uniform translation s that results in overlaps. 
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FIG. 7. Geometrical model construction in Conpan: (a) 
Two segments of the array, (b) First-layer sections generated 
by Conpan for the two segments, (c) Second-layer island sec- 
tions (shaded), partially overlapping first-layer sections, (d) 
Shape of first-layer islands, (e) Shape of second-layer islands 
overlapping first-layer, (f) First- and second-layer islands 
paneled for capacitance calculation. 



1. Overlap Detection 

Since a single second-layer section can overlap more 
than one first-layer section, Conpan detects the overlaps 
using a recursive detection algorithm. To begin, each 
second layer section is compared against each first layer 
section to detect overlaps. When an overlap is found, 
the second layer section spawns two daughter sections, 
one overlapping and one not. The axis and coordinate of 
the split are stored in the the mother section, along with 
pointers to the daughter sections. The mother section 
becomes a placeholder, used only to keep track of the 
relationship among its daughter sections. 

The non-overlapping daughter is then is compared 
against the remaining first-layer sections to find other 
overlaps. If there are more overlaps, the daughter spawns 
a pair of sub-daughters, and so on. The recursive process 
stops when no new overlaps are found. The daughter sec- 
tions that remain undivided are called "final daughters" . 
Figure || gives a schematic view of the recursive overlap 
detection process for a single island. 
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FIG. 8. Recursive overlap detection in Conpan for a circuit 
fragment where a second-layer island S overlaps two first-layer 
islands F 1 and F 2 . (a) Top view, with overlap areas shaded, 
(b) Comparing S to all first-layer islands, overlap of F 1 is 
detected. S is split into two daughter sections, S a and St, at 
y = y* . (c) Comparing St, to all first-layer islands, overlap 
of F 2 is detected. St is split into Sba and Sbb a t V = y%- 
(d) Daughter sections S a , Sba, and Sbb, having no further 
overlaps, are used to build the 3D structure shown in (e). 

In addition to dividing up the second-layer sections 
to account for overlaps, Conpan also splits first-layer 
sections along the line where they are overlapped (see 
Fig. 0f ) . This allows the edges of panels facing each other 
across a junction to line up, facilitating convergence in 
capacitance calculations. 



2. 3D Representation 

Once all the overlaps have been found, Conpan can 
begin to create the 3-D model of the circuit. Each fi- 
nal daughter section becomes the base of a "block", a 
rectilinear solid representing part of a conductor. The 
heights of the two layers are specified by the parame- 
ters h\ and hi. First-layer blocks and non-overlapping 
second-layer blocks have their base at z = 0. Overlap- 
ping second-layer blocks have their base at z = hi + 1, 
where t is the thickness of the gap between overlapping 
islands that represents the tunnel junction. 
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Finally, each block surface is divided into panels. The 
goal is to divide the surfaces in such a way that an accept- 
ably accurate capacitance calculation can be performed, 
within the limits of available computer memory and cal- 
culation time. The division process is guided by an in- 
put parameter a, the goal panel length. The surface of a 
block with length Li along axis i is divided into the num- 
ber m of divisions that brings Lj/n, closest to a. Once 
the block surface has been divided along both its axes, 
the resulting grid of panels is written to a panel file for 
input to FastCap. Each panel is stored simply as a quar- 
tet of x, y, z coordinates, one for each of the four corners, 
together with the number of the node it belongs to. 



B. Conpan Input Parameters 

1. Junction Thickness 

In the physical circuit, the tunnel barriers separating 
the islands consist of AlO^, with unknown x and thick- 
ness tj. From literature data on similar junctions [ p^[ , 
we expect a dielectric constant ej ~ 4 and tj ~ 1 nm. 
FastCap can handle dielectric surfaces much as it han- 
dles conductors - by dividing them into panels. However, 
each additional dielectric panel demands more computer 
memory and calculation time. Since tj is much smaller 
than the transversal dimensions in all junctions, the elec- 
tric field configuration outside the junctions does not de- 
pend strongly on their internal geometry. Therefore, we 
avoided modeling the junction dielectrics explicitly by re- 
placing them with uniform free-space gaps (e = 1) with 
the effective thickness t = tj/ tj . This effective thickness 
was adjusted to make the junction specific capacitance 
match the standard experimental value 4.5/iF/cm 2 typi- 
cal for the Al/AlOz/Al junctions with tunnel conductiv- 
ity in our range (~ 10 5 S/cm 2 ) |l5 16|. 



2. Line Widths 

The effective line width w of islands (and of the nar- 
row parts of externals) is difficult to measure directly, 
because of its small magnitude (see Fig. ||b). We de- 
termined w by requiring that the simulated inverse self- 
capacitance of the electrometer island (C 1 ^ 1 14 ) match its 
experimentally measured value. We derive C 1 ^ 1 14 from 
the maximum value of the electrometer Coulomb block- 
ade threshold voltage Ut, as seen in electrometer I-V plots 
(Fig. |): 
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FIG. 9. Current-voltage traces for the electrometer, 
(SiCh/Si substrate, T=35mK, superconductivity suppressed 
by magnetic field), current biased, taken while rapidly vary- 
ing Vs . The Coulomb blockade voltage Ut varies with Vs from 
to a maximum value of eCj 4 ^ . 

Ci^m (conventionally known as Cs) in turn depends 
on w because it is is dominated by the electrome- 
ter junction capacitances, which increase monotonically 
with w. The experimentally measured values for C 14 1 14 
were (4.6 x 10 _16 F) _1 for the circuit on Si (sample 
#LJS011494B) and (2.7 x KT^F)- 1 for the circuit on 
Si0 2 /Si substrate (sample #LJS01I494A). The island 
width w used in simulation, as determined from Eq. (|^) , 
was 30 nm for Si and 42 nm for Si02/Si. Both of these 
values are consistent with the values expected from fab- 
rication parameters and from AFM and SEM imaging of 
the samples. 

W , the width of the wide parts of the externals, is 
specified as 1 /jm in our mask files, and can be accepted 
at "face value" because it is large compared to the scale 
of geometrical uncertainty in the circuit, and because the 
wide parts of the externals are all far (several /im) from 
the islands. 



3. Layer Heights 

The heights of the two layers (hi and I12) are deter- 
mined with a quartz monitor in the deposition unit dur- 
ing fabrication. In our case, these heights were measured 
to be 30 and 50 nm (± 10%), respectively. 



C. Substrate 



To calculate the effects of the substrate on circuit ca- 
pacitances, FastCap requires a paneling of the complc- 
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mentary image of the "footprint" of the nodes, because 
panels representing the dielectric/metal interface (parts 
of the substrate covered by nodes) must be treated differ- 
ently than panels representing the dielectric/air interface 
(the exposed substrate). In a manner analogous to that 
used for conductor panels (see Section IV), one could in- 
vestigate various methods of paneling the complementary 
substrate image in order to minimize the number of pan- 
els, while yet retaining an acceptable level of accuracy. 
Such a paneling algorithm itself is not simple to create. 

We avoided this problem through an old calculational 
trick in electrostatics - the image method. A modified 
version of FastCap was created, called ImageCap, which 
can simulate the effects of a single- or double-layer sub- 
strate by creating a set (or multiple sets, in the double- 
layer case) of image panels. 
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FIG. 10. Image panels in ImageCap: (a) single substrate, 
(b) double substrate. 

In the single substrate case, each image panel is formed 
by reflecting the original panel about the plane of the sur- 



face of the substrate (Fig. |lOja) . For the purposes of cal- 
culating the electrostatic potential above the substrate, 
the charge on the image panel is 
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1 



(G) 



where q is the charge on the original panel, and e is the 
relative dielectric constant of the substrate. 

For a substrate covered by an oxide of thickness H, 
an infinite series of image charges is required for an ex- 
act representation of the electrostatic effect of the sub- 
strate (Fig. |To|b). However, the distance from the original 
charge to each successive image charge increases linearly, 
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while the value of each successive image charge decreases 
exponentially, 



92,,: = 4/3 



i = 1,2. 
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Here t\ and e 2 are the dielectric constants of the surface 
oxide layer and the bulk substrate, respectively. For our 
circuits, we have accepted the table values e = 12.1 for 
the bare Si substrate and ei = 4.5, e 2 = 12.1 for the 
Si0 2 /Si substrate (a = 0.64, /3 = -0.46). The resulting 
expression for the double-layer image charges, 



Q2,i 



-0.36 x (0.29) 



Of-i), 



i = 1,2. 



(10) 



shows that qiA is already down by three orders of magni- 
tude from the original charge. In our calculations, adding 
image levels beyond g 2j 4 made no difference to the result, 
within a relative error (of the largest self-capacitances) 
below - 10" 4 . 



IV. CAPACITANCE MATRICES 

A. Matrix Structure 

Using the circuit panels generated by Conpan, Image- 
Cap generates the capacitance matrix for the circuit. Im- 
ageCap adds the effect of image panels when calculating 
potentials, and uses no multipole acceleration; otherwise, 
its algorithms are the same as in FastCap @. First, the 
inverse capacitance matrix for panels is calculated and 
inverted. Each element dj in the capacitance matrix for 
nodes is then formed by summing all the panel capaci- 
tance matrix elements linking nodes i and j. The charges 
and potentials on the nodes are related by 
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so that Cij is numerically equal to the amount of charge 
induced on node j when node i is held at unit potential 
and all other nodes have zero potential. 

C is an N x N matrix, where TV = N e +Ni, and N e and 
Ni are the numbers of external nodes and island nodes in 
the circuit, respectively. Ordering all the external nodes 
before the island nodes, we can write C in terms of sub- 
matrices: 



C = 



<g> 


-c 


-c T 


c 
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Here C is the symmetric Ni x N matrix of island- 
island capacitances and C is the N e x Ni matrix of 
external-island capacitances (with elements defined pos- 
itive, by convention). External-external capacitances 
(represented above by the ®) are not needed for our 
simulations. 

The matrices calculated for our circuits are shown in 
Tables I and II. Note the up/down alternation of mu- 
tual capacitances along the array for the circuit on Si - 
e.g., C\,i, the capacitances linking external node 1 to the 
islands. For example, although island 7 is closer to ex- 
ternal 1 than is island 8, C x ,i ^ 0.030 x KT 16 F is smaller 
than Ci,s — 0.045 x 10~ 16 F (similarly for islands 9 and 
10). In C, we see that || C 6>9 1| is smaller than || Cs,xo || ? etc. 
This phenomenon reflects the influence of the silicon sub- 
strate, which, due to its high dielectric constant (e ~ 12), 
links externals to the first-layer islands (which lie flat on 
the substrate) more strongly than to the second-layer is- 
lands (which lie partly on top of the first-layer islands). 
The capacitances for the circuit on Si02/Si do not show 
these oscillations as strongly, as we would expect from 
the smaller permittivity (~ 4.5) of SiC>2- 



B. Model Accuracy 

Our model contains three main simplifications related 
to computational constraints, each of which introduces 
error into our capacitance matrix calculations. 



1. Free-space Junctions 

As noted above, we calculate C using free-space junc- 
tions of thickness t instead of dielectric junctions of thick- 
ness tj. (Although initially this approximation was in- 
tended for convenience, it later became a necessity as Im- 
ageCap does not handle explicit dielectric panels.) The 
error involved in this approximation was estimated by us- 
ing FastCap to model a chain of islands in two ways: with 
explicit dielectric junctions and with free-space junctions. 
Results for an 8-island chain, with effective dielectric 
thickness chosen to make the island self-capacitances in 
both models the same, indicate that the error involved 




in this approximation is below 1% for junction-linked is- 
lands, and between 1% and 4% for non-junction-linkcd 
islands. 



2. Paneling 

In calculating capacitances, FastCap/ImageCap as- 
signs a uniform charge distribution to each panel. Hence, 
its accuracy depends on how well the paneling follows 
changes in change distribution on the node surfaces. 
Clearly, the denser the paneling, the better the repre- 
sentation of changes in charge distribution. However, 
panel density is effectively limited by available computer 
memory. For example, a FastCap simulation with 5000 
panels typically requires more than 128 MB. ImageCap 
uses even more memory, since it calculates all panel in- 
teractions directly. We investigated the dependence of 
calculated capacitance on paneling density for a simple 
two cube system (Fig. [TT|a). 
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FIG. 11. Two cube system with L/t = 10: (a) Panelings 
with 1, 3, and 9 panels/side, (b) Capacitances as a function 
of panel density (q — C(j>)- 

The results (Fig. [ljjb) suggest that a non-uniform 3x3 
grid (with a 1/10 ratio of edge panel length to central 
panel length, reflecting the peak in surface charge near 
the edges) for the smaller, roughly square-shaped node 
faces (Fig. Qf) is sufficient to calculate capacitances with 
an error below 10%. This is essentially how we paneled 
roughly square-shaped island surfaces. For longer faces 
(Fig. |l2|a) we used a larger number of divisions along 
their length. In an islands-only test circuit, increasing 
the total number of panels from ~ 2000 (corresponding 
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to the 3x3 grid for roughly square-shaped surfaces) to 
6000 resulted in less than 1% changes in island-island 
capacitances. Thus we believe that the total error in 
island-island capacitances due to finite panel density is 
perhaps only ~ 1%. 




FIG. 12. Paneling used in our calculations: (a) Trap islands 
(nodes 12 and 13). (b) Thin parts of externals, (c) Wide parts 
of externals. 

To reduce the number of panels in the model, the two 
layers of an external are fused into one where they over- 
lap. The error involved in this simplification is negligible. 
In addition, the narrow parts of the externals were di- 
vided along their length without edge panels (Fig. |l2|b). 
This simplification was found to cause an error in island- 
external capacitance of ~ 5% when the island and the 
external are connected by a junction (Fig. |l3|), and ~ 1% 
otherwise. Finally, wide parts of the external leads 
were represented by only their top and bottom surfaces 



(Fig. [T^c), again to save panels. Since the width to height 
ratio W/h ~ 12, the error introduced by this simplifica- 
tion is negligible. The top and bottom surfaces are di- 
vided according to the 3x3 type scheme described above 
for islands. Despite the large size of the resulting panels, 
the error involved in this simple paneling appears to be 
~ 1%. 
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FIG. 13. Effect of panel density on the calculated is- 
land-external capacitance. The test circuit had only islands 
14 and 15 and narrow parts of externals. 

3. Lengths of Externals 

The calculated capacitance values depend on the 
lengths of the external wires used in the model. In gen- 
eral, island-external capacitances increase with external 
length, at the expense of island stray capacitance (ca- 
pacitance to a ground at infinity); the self-capacitance 
of islands does not change appreciably. To measure the 
error introduced by cutting off the externals at a given 
length, we have calculated capacitance matrices for test 
circuits with varying external lengths (Fig. [l4]). These 
circuits consisted of only one island and only the wide 
parts of the five externals. As a result, the error induced 
by cutting off externals in these test circuits should be 
proportionately larger than the error in the complete cir- 
cuits. Still, the test circuits indicate that the error in- 
volved in cutting of the circuit at a radius of 20 /im (as 
in our final versions of the complete circuits) was less 
than 2%. 
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FIG. 14. Effect of lead length cutoff on the calculated is- 
land-external capacitance. The test circuit had only island 12 
and wide parts of externals. 



4- Total Error 

Considering the error caused by the above simplifi- 
cations in the calculation of C itself, it seems safe to 
say that the the combined error for any given calculated 
capacitance matrix element was less than 10%. Note 
that we are not yet considering how well the geomet- 
rical model corresponds to the physical circuit (see Sec- 
tion jv|. 



V. SIMULATED AND EXPERIMENTAL 
RESULTS 

We have calculated most properties of our circuits us- 
ing MOSES, the single-electron circuit simulation program 
|p7| . This program uses a Monte Carlo algorithm to sim- 
ulate arbitrary SET circuits within the framework of the 
orthodox theory of single-electron tunneling |l],|| . MOSES 
needs to know the capacitance sub-matrices C and C and 
the conductances of all tunnel junctions. The resistance 
of two electrometer junctions connected in series can be 
extracted from the slope of the experimental dc I — V 
curve of the electrometer at high voltage (V ^> e/Cs). 
From this measurement, we calculated tunnel conduc- 
tance per unit area. Conductances of all other junctions 
in the circuit were then calculated by assuming that their 
conductance is proportional to their nominal area. This 
assumption may only be accurate to an order of magni- 
tude; however, most of the results discussed below per- 
tain to stationary properties of the system, and are thus 
unaffected by deviations in conductance. 



A. General Electrostatic Relations 

Solving the matrix equation ( [TT| ) for the island poten- 
tials (pi, with our definition ([l2|) of the capacitance matrix 
we get 

& = E C ^9i + Qj) > Qi = £ Ck i Vk > ( 13 ) 

j£isl kGext 

or, in a different form, 
4>i = E C ij 1( b + aikVk ' aik = E Cij^kj , 

(14) 

where V% are the external potentials. These relations al- 
low us to establish useful relations between changes in 
the external potentials {Vfc} and the charge state of the 
islands {qi}, and the dynamics of the system as deter- 
mined by the island potentials 



B. Electrometer 

Let us apply these relations, in particular, to the is- 
land of the single-electron transistor (number 14 in our 
notation, see Fig. |l|) serving as the electrometer. Exper- 
imentally, we measure the dc voltage U3-4 between the 
"source" and "drain" of the transistor (externals 3 and 
4) under a small (~ 100 pA) dc current bias. If the tem- 
perature is small enough (AsbT <C e 2 C 1 ^ 1 14 ), the voltage 
U3-4 in such an experiment closely follows the threshold 
Ut of the Coulomb blockade of the transistor - see Fig. ||. 

It is well known (see, e.g., Refs. that the thresh- 

old is determined by the effective background charge Q a 
of the transistor island, which may be defined as 

014|£/ 3 =£/ 4 =O = ^14^14 (914 + Qo) ■ (15) 

Comparing ( |l5| ) and ( jl4| ) above, we obtain in our nota- 
tion 



1 



Qo — 1 



J'#14 

E G rlj1i + E aikVk 

jGisl k£ext 



(16) 



Eq. ( |l6] ) allows us to find the theoretically expected 
variation of Q a due to any changes in the system. On 
the other hand, the threshold voltage is an e-periodic 
function of Q D , and its maximum amplitude is expressed 
by Eq. (|^) (for the case when the two transistor junction 
capacitances are the same). Thus, after we measure the 
experimental value of (Ui) m ax, we can express the change 
in the effective charge Q a via the observed variation in 
U,: 

^\^t Jmax 

We have applied this approach to compare experiment 
and theory for two samples (#LJS011494A with Si0 2 /Si 
substrate and #LJS011494B with Si substrate). 
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C. High-T electrometer response 



D. Trap phase diagram 



We can readily measure AVi (i = 1,2,5), the change in 
external voltage corresponding to one period of the oscil- 
lating threshold voltage (Fig. At k B T > 0.1e 2 /C E 
(experimentally, T > 0.5K), thermal activation of elec- 
trons smears the Coulomb blockade effects and makes the 
junctions essentially transparent to tunneling, while the 
periodic response of the electrometer is still visible up to 
k B T ~ 0.3e 2 /C s (T ~ 1.5K). Thus, the measured val- 
ues of AVi depend only on the circuit geometry and are 
essentially independent of the properties of the junctions. 
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FIG. 15. A typical experimental dependence of electrome- 
ter dc voltage U3-4 on gate voltage V5 for a circuit comprising 
only an electrometer. 

MOSES is not a useful tool for directly modeling high-T 
behavior, as the number of jumps involved would be ex- 
tremely high. However, we can simulate high-T behav- 
ior in MOSES by specifying very high external voltages 
Vi while keeping temperature low (say, T = 0). Un- 
der these high voltage conditions, the islands are flooded 
with extra electrons, and the tunnel junctions become 
effectively transparent to tunneling, just as in the high 
temperature case. Thus we simply apply an external volt- 
age Vi S> e/Cs, measure how many electrons enter the 
electrometer island, and find the ratio of voltage Vi to 
electrometer charge qu. 

Table III shows values of the ratio AVi for simulated 
circuits and for experimental circuits averaged over sev- 
eral nominally identical samples. For the experimental 
values, the uncertainties given reflect the spread of among 
the samples. For the simulated values, the uncertainties 
given reflect the ~ 10% error in calculated values, as de- 
scribed in Section IV. The simulated values are all lower 



than the experimental ones (with the exception of AV5 
on Si02/Si), differing by as much as 50%. The agreement 
is somewhat better for the circuits on SiC^/Si. 



The simplest measurable characteristic involving 
single-electron charging of the trap is its phase diagram 
(Fig. |l6|), which reflects changes in the charge states of 
the array and trap as a function of the drive voltage V\ . 
In MOSES, we can directly view the charge state of each 
island in the array and trap as we vary V\ , as well the re- 
sulting change in Q a . In the physical circuit, however, we 
can only measure the response of the electrometer 

and reduce it to the changes in Q using Eq. [l7| 
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FIG. 16. Electrometer phase diagram: Qo as a function of 
trap drive voltage. V2=8.2mV, sample on Si02/Si. Crosstalk 
from drive voltage to electrometer has been subtracted, leav- 
ing only influence of trap charge state. In experiment, su- 
perconductivity in aluminum is suppressed by a 2T magnetic 
field. 

Figure ^ shows experimental and simulated electrom- 
eter response to ramping the trap drive voltage V\ up and 
down over a period of several minutes. In both cases, 
the effects of the crosstalk between external 1 and the 
electrometer have been removed. In the experiment, the 
crosstalk is cancelled by feeding the electrometer gate 
(node 5) with a voltage V5 = —aVi, with the coefficient 
a adjusted to make the phase diagram plateaus horizon- 



ff 



tal. In simulations, MOSES accomplishes the same effect 
by subtracting A^i4 = 01,14 Vi from the electrometer is- 
land potential. 

Horizontal plateaus in Fig. ^ correspond to particular 
charge states of the system (trap + array), while verti- 
cal jumps correspond to changes of charge state. Thus, 
the hysteretic loops are regions of bi/multi-stability. The 
blow-up of the theoretical curve (Fig. |lj]) indicates the 
states for several plateaus. In particular, notice that the 
largest plateaus correspond to states that are most sta- 
ble because the array is either charged uniformly (one 
electron on each island, for example) or in a regular al- 
ternating pattern such as 1-0-1-0 (Fig. |l8|). The smaller 
plateaus correspond to more complex charge states which 
are less stable. 
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FIG. 17. Closeup view of experimental and simulated 
phase loops. Charge vectors for islands 6-13 on various 
plateaus in simulated phase diagram are indicated in brackets. 
In simulation, q = was assumed. 
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FIG. 18. Simulated phase loops at a higher range in Vi, 
showing wider plateaus around [1010123]. 

The experimental phase diagram bears a qualitative re- 
semblance to the theoretical one, with somewhat shorter 
plateaus, though the order of magnitude is the same 
(~ 2 mV for major plateaus). Simulated phase diagrams 
with randomly selected qo show shorter plateaus than the 
qo = phase diagram (see Secjv| below). 

In Fig. [H], the large jumps in Q a correspond to a single 
electron entering the trap: 5Q = (AQ c ,)| e _ >tr . Using 
Eq. (plf), we can also express the simulated value of SQ 
as 



SQ S o im = 7 ^ ± e, 

C14 14 



(18) 



where tr — 12 or 13, depending on which trap island 
the electron stops in. For comparison with experimental 
results, we take the average of the two possible values. 
The results are shown in Table IV. The difference be- 
tween simulated and experimental values for Si is within 
the estimated geometric calculation error (10%), while 
the value for SiC^/Si is not. 



E. Plateau dependence on V2 

For a given plateau, the switching voltages V\ = V± 
depend on the "ground" voltage V2 (see Fig. |l|). In the 
simplest model, with no stray capacitances, (see, e.g., 
Ref. ||) the charge state of the system depends only on 
the voltage V = V\ — V%. In that model, the depen- 
dences V±{V2), corresponding to changes in the charge 
state, would form parallel 45° lines in the [V\, V2] plane. 
In reality, however, stray capacitances of the islands to 
"infinity" (i.e. to a distant common ground) make the av- 
erage potential (V\ +V2V2 of the system relevant as well. 
As a result, the region corresponding to each charge state 
acquires a shape similar to a stretched diamond (Fig. fl9l ). 
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FIG. 19. Experimental and simulated threshold voltages 
V± as functions of V2: Si substrate, superconductivity sup- 
pressed. 

Simulations using MOSES show that the diamond shape 
results from the alternation of two types of electron trans- 
port that switch the charge state. At the I0W-V2 end of 
the diamond, the charge state switches with the transfer 
of an electron in/out of the trap (see the energy diagram 
in Fig. ||a). However, at the high- V2 side, the barrier 
for holes to enter or exit is lower than for electrons (cf. 
Fig. ||b). Near the sharp ends of the diamond, the crit- 
ical transport may be even more complex (e.g., creation 
of an electron-hole pair inside the array, with the sequen- 
tial motion of its components apart, one into the trap, 
and another into the external electrode). In these re- 
gions, however, the plateau corresponding to the charge 
state of the trap is already small and virtually disappears 
among numerous plateaus corresponding to various in- 
ternal charge states of the array (Fig. |l7]). Figure [l^ 
shows that while the diamond shape of the charge state 
in [VijVy is well reproduced in experiment, the simu- 
lated width || V+ — V—\\ of the bistability region in V\ is 
roughly twice the experimental value. 

For each V2, there is one value of V\, called V eq , at 



which the energy barrier is the same for an electron to 
tunnel into or out of the trap fL8fl . A good measure of 
the relative influence of the two external voltages on the 
trap is the derivative 



eg 



«2,tr 



dV 2 



1 - ax.tr 



(19) 



where tr — 12 or 13, depending on which trap island ac- 
tually traps the electron for a given (Vx,V2)- The two 
values are typically within 5% of each other, and we take 
their average when comparing simulated and experimen- 
tal results. In the experimental data, we define the av- 
erage V eq by bisecting the diamond shape in the graph 
(Fig. |l^). dVeq/dVz is essentially a geometric property 
of the circuit and should not depend on thermal activa- 
tion or cotunneling. As Fig. [l9] shows, the simulated and 
experimental values are very close. 



F. Energy barrier 

At Vi = (V2 ) , we can measure the energy barrier 
AW experimentally by measuring trapping lifetime as 
a function of temperature (for experimental details, see 
Ref. [0). The Arrhenius law for lifetimes gives 



r L oc exp(AW/kT), 



(20) 



so that plotting log(rt) vs. 1/T gives us AW. Dynami- 
cal simulations |lq] have shown that ( po|) is virtually un- 
affected by cotunneling for relatively high temperatures 
(~ 100 mK and above). In simulation, MOSES allows 
us to measure AW directly. Figure ^ shows the depen- 
dence of the trap energy barrier AW on the bias voltage 
V2, for the circuit on the SiC^/Si substrate. The sim- 
ulated energy barrier profile peaks at roughly the same 
value of V2 as in the experiment, and the peak barrier 
value is within ~ 10% of the experimental value. How- 
ever, the simulated peak is sharper. 
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FIG. 20. Experimental and simulated dependence of bar- 
rier height on Vi (Si02/Si substrate). 



VI. DISCUSSION 

Let us first discuss the results independent of the 
single-electron charging effects: the transistor response 
5Q to a single electron entering the trap, the oscillation 
periods AVi, and the slope dV eq /dV-2- The differences be- 
tween simulated and experimental values for SQ are 6% 
and 23% for the Si and SiO^/Si substrates, respectively. 
Values for AVi do not agree as well: differences between 
simulated and experimental values range from 15 to 38% 
for the trap on SiC^/Si, and from 29 to 50% for the trap 
on Si. We had experimental data for dV eq /dV2 only on 
Si. Here the difference between experimental and simu- 
lated results was ~ 12%. These numbers suggest how 
well our geometrical model corresponds to the physical 
circuit (the accuracy of the orthodox theory and of the 
MOSES simulator is presumably much higher). 

The most obvious idealization involved in our geomet- 
ric modeling is that the islands created by Conpan are 
rectilinear and uniform. Even at the limited resolution 
of an AFM image (Fig. |]b), the contours of the fabri- 
cated circuits appear rounded and irregular on a scale of 
~ 10 nm. This is to be expected, due to the relatively 
large grain size of evaporated Al (~ 50 nm, compara- 
ble to the line width w) and the stochastic nature of the 
grain growth process. Most capacitance matrix elements 
should not depend strongly on small details of the island 
shape. However, irregularities in the shape of overlapping 
islands may change the area, and thus the capacitance, 
of the junctions linking them. 

All other results involve single-electron charging ef- 
fects. Here the difference between the theory and experi- 
ment is larger - typically by a factor of 2, and sometimes 
larger. We believe that the most important origin of this 
difference is the set of background charges q . The Si sub- 
strate is capable of trapping charged impurities near the 
circuit islands. The result of these impurities is that the 
charge on island i effectively changes from q\ to q\ + qoi . 
These charges may furthermore be capable of thermal 
migration over time. 

Simulated plots of the electrometer response to trap 
charging with three randomly selected <fo are shown in 
Fig. EU It appears that the wide (~ 4 mV) steps near 
Vx = in the qo — plot are not stable to variations in 
go: in most plots with random qo, as in the experimental 
plot, all step widths are less than 3 mV. 
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FIG. 21. Simulated response of electrometer to ramps of Vi 
for go = (solid line) and three randomly selected background 
charge vectors qo. Plots are shifted vertically because back- 
ground charge on the electrometer island (node 14) effectively 
shifts Qo- 

To summarize: we have developed an automated way 
to construct simplified models of experimental single- 
electron devices and circuits with metallic islands, and 
we have compared the properties of model single-electron 
traps with those of real traps. The observed differ- 
ences between simulation and experiment may be at- 
tributed to random deviations of the physical struc- 
tures from their nominal size and shape, and to ran- 
dom background charges created by charged impuri- 
ties. Future work of interest may include more precisely 
defined single-electron devices, using better fabrication 
technology, and the extension of quantitative modeling to 
semiconductor-based single-electronic circuits and hybrid 
single-electronic / conventional transistor logic circuits. 
As single-electronics evolves into a mature technology, 
such modeling will be essential. 
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3.7254 




















7 


-1.7587 


3.6101 
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-0.0223 


-1.7587 


3.7259 
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-0.0038 


-0.0084 


-1.7587 


3.6103 














10 


-0.0056 


-0.0038 


-0.0222 


-1.7586 


3.7261 












11 


-0.0016 


-0.0012 


-0.0037 


-0.0084 


-1.7585 


3.6104 
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-0.0062 


-0.0036 


-0.0117 


-0.0075 


-0.0352 


-1.7678 


5.7693 








13 


-0.0032 


-0.0019 


-0.0052 


-0.0033 


-0.0096 


-0.0118 


-3.4910 


3.9670 






14 


-0.0020 


-0.0010 


-0.0027 


-0.0014 


-0.0034 


-0.0016 


-0.0312 


-0.0228 


4.6482 




15 


-0.0018 


-0.0009 


-0.0024 


-0.0012 


-0.0031 


-0.0015 


-0.0277 


-0.0349 


-0.2780 


0.9550 



TABLE I. Capacitance matrix generated by ImageCap for circuit on Si substrate. Rows 1-5 belong to C, rows 6-15 belong 
to C. By convention, all CV, > 0. All values in 10~ 16 F. 
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7 


-1.2370 


2.5324 


















8 


-0.0070 


-1.2370 


2.5501 
















9 


-0.0017 


-0.0067 


-1.2370 


2.5321 














10 


-0.0015 


-0.0017 


-0.0070 


-1.2369 


2.5500 












11 


-0.0006 


-0.0009 


-0.0017 


-0.0066 


-1.2369 


2.5323 










12 


-0.0016 


-0.0015 


-0.0033 


-0.0035 


-0.0111 


-1.2418 


3.7851 








13 


-0.0009 


-0.0010 


-0.0016 


-0.0019 


-0.0034 


-0.0084 


-2.3386 


2.5605 






14 


-0.0005 


-0.0004 


-0.0007 


-0.0006 


-0.0009 


-0.0007 


-0.0112 


-0.0087 


2.6644 




15 


-0.0005 


-0.0004 


-0.0007 


-0.0006 


-0.0009 


-0.0007 


-0.0100 


-0.0145 


-0.1097 


0.4104 



TABLE II. Capacitance matrix generated by ImageCap for circuit on SiCh/Si substrate. Rows 1-5 belong to C, rows 6-15 
belong to C. All values in 10~ 16 F. 



1G 



Node i 


A yexp 


Si substrate 

A ysim 


Difference 


A yexp 


Si02/Si substrate 

Avr m 


Difference 


1 


16(2) 


8(1) 


50% 


32(5) 


20(2) 


38% 


2 


24(1) 


17(2) 


29% 


54(2) 


45(5) 


15% 


5 


14(2) 


7(1) 


50% 


38(6) 


51(5) 


34% 



TABLE III. Values for AVi in mV. 



Substrate 


5QT P 


SQV m 


Difference 


Si 


0.064e 


0.060e 


6% 


Si0 2 /Si 


0.059e 


0.045e 


23% 



TABLE IV. Response of electrometer to single electron entering the trap 
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